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Selma Koghee^'^, Lih-King Lim^, M.O. Goerbig^, and C. Morais Smith^ 
^Institute for Theoretical Physics, Utrecht University, 
Leuvenlaan 4, 3584 CE Utrecht, The Netherlands and 
^ Laboratoire de Physique des Solides, CNRS UMR 8502, Universite Paris-Sud, 91405 Orsay, France 

(Dated: February 1, 2012) 

Inspired by the recent creation of the honeycomb optical lattice and the realization of the Mott 
insulating state in a square lattice by shaking, we study here the shaken honeycomb optical lattice. 
For a periodic shaking of the lattice, a Floquet theory may be applied to derive a time-independent 
Hamiltonian. In this effective description, the hopping parameters are renormalized by a Bessel 
function, which depends on the shaking direction, amplitude and frequency. Consequently, the hop- 
ping parameters can vanish and even change sign, in an anisotropic manner, thus yielding different 
band structures. Here, we study the merging and the alignment of Dirac points and dimensional 
crossovers from the two dimensional system to one dimensional chains and zero dimensional dimers. 
We also consider next-nearest-neighbor hopping, which breaks the particle-hole symmetry and leads 
to a metallic phase when it becomes dominant over the nearest-neighbor hopping. Furthermore, 
we include weak repulsive on-site interactions and find the density profiles for different values of 
the hopping parameters and interactions, both in a homogeneous system and in the presence of 
a trapping potential. Our results may be experimentally observed by using momentum-resolved 
Raman spectroscopy. 



INTRODUCTION 



The study of Dirac points, i.e. the contact points be- 
tween different energy bands with an approximate linear 
dispersion relation, has become a major issue since the 
experimental breakthrough in graphene-based electron- 
ics [1][2]. Indeed, the low-energy electronic properties of 
graphene are governed by a pseudo-relativistic 2D Dirac 
equation for massless fermions situated at the K and K' 
corners of the Brillouin zone [3 . The Dirac points are 
topologically protected and a gap is opened only when 
the inversion symmetry of the lattice or the time-reversal 
symmetry are broken. 

The possibility to generate topological phase transi- 
tions in graphene-like systems has recently attracted a 
great deal of attention. Within a tight-binding descrip- 
tion, an anisotropy in the nearest-neighbor hopping pa- 
rameters makes the Dirac points move away from the 
high-symmetry K and points and, under appropri- 
ate conditions, merge at time-reversal invariant points in 
the first Brillouin zone [4]-[6]. Most saliently, this merg- 
ing of Dirac points is associated with a topological phase 
transition between a semimetallic phase and a gapped 
band-insulating phase. An experimental investigation of 
the merging transition in graphene turns out to be prob- 
lematic, since in order to appropriately modify the hop- 
ping parameters, an unphysically large strain needs to be 
applied to the graphene sheet [7 . 

An alternative system for the study of such topolog- 
ical transitions is that of ultracold atoms trapped in a 
honeycomb optical lattice. Since the seminal realization 
of the superfluid-Mott-insulator transition in the Bose- 
Hubbard model, ultracold atoms in optical lattices have 
become promising systems to emulate condensed-matter 
physics. Indeed, the lattice geometry, the dimensionality. 



the atomic species, as well as the interactions can be engi- 
neered with a high degree of precision. The more involved 
triangular and honeycomb geometries were recently real- 
ized experimentally and exotic correlated states of mat- 
ter have been observed experimentally [8 or predicted 
theoretically [MS. 

The application of a time-periodic perturbation on the 
optical lattice introduces yet another parameter scale 
into the system. A periodic shaking of the optical lattice, 
up to the kHz frequency range, has been implemented by 
placing one of the mirrors used to create the optical lat- 
tice on a piezoelectric material, such that the mirror can 
be moved back and forth in the direction of the beam 
[131 • The Floquet formalism shows that the hopping 
energy of the atoms in the shaken lattice is renormalized 
by a Bessel function, as a function of the shaking fre- 
quency and amplitude, thus allowing both the magnitude 
and the sign of the hopping parameter to change. This 
rather counter-intuitive phenomenon, as compared to the 
standard tight-binding physics, has been experimentally 
observed in a one-dimensional cold-atomic system 

In this paper, we consider ultracold fermions trapped 
in a shaken honeycomb optical lattice. Within the Flo- 
quet formalism, we derive an effective Hamiltonian that 
generalizes that of a graphene-like material under strain. 
In particular, we find that the alignment and merging 
of Dirac points in momentum space are now accessible 
with ultracold fermions in the shaken optical lattice and 
the phase diagram consists of various phases of the corre- 
sponding solid-state system that are otherwise difficult to 
realize. Furthermore, by taking into account a Hubbard- 
like interaction for spinful fermions, we study the density 
profiles for the homogenous and the trapped gas within 
a Hartree-Fock theory. 



The outline of this paper is the following: in Sec. [HA 
we introduce the time-dependent Hamiltonian and in 
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Sec.|IIB| we derive the time-independent one, by apply- 



ing the Floquet formahsm. In Sec. Ill we investigate the 
merging and ahgnment of Dirac points, when the optical 
lattice is shaken along specific directions. The descrip- 
tion is extended to include interactions in Sec.[TVl where 
we derive the dependence of the density on the chemical 
potential. Implications of our results for experiments are 
discussed in Sec.jVl Finally, our conclusions are presented 
in Sec.lVD 



II. THE SHAKEN HONEYCOMB LATTICE 



energy gained in hopping to the nn and next-nearest- 
neighbor (nnn) sites, respectively, and /i is the on-site 
energy. We remark that the nnn hopping is taken into 
account because the nn hopping may be rendered vanish- 
ingly small in the effective time-independent description. 
In this regime, the nnn may become the dominant kinetic 
term. In a square lattice, where the potential is separable 
in independent Cx and components, the nnn hopping 
is identically zero [15]. However, the nnn hopping can be 
nonzero in the honeycomb lattice, since its potential is 
not separable in Cx and Cy components. Nevertheless, it 
may be expressed as the sum of two triangular lattices. 



In this section, we derive a time-independent effec- 
tive description for ultracold atoms trapped in a peri- 
odically shaken honeycomb optical lattice by utilizing a 
Floquet theory. For simplicity, we focus on a system 
of single-component fermionic atoms and consider only 
single-particle terms in this section. The results from 
the Floquet theory are valid for fermionic atoms with in- 
ternal degrees of freedom as well as for bosonic atoms. In 
particular, the hyperfine state of fermionic atoms, play- 
ing the role of an effective spin- 1/2 degree of freedom for 
electrons, will be considered when interaction effects are 
taken into account in Sec. IV. 



A. Time-dependent Hamiltonian 

In the tight-binding limit, the system of ultracold 
fermionic atoms trapped in a 2D shaken honeycomb op- 
tical lattice can be described by the Hamiltonian 

H{t) = Ho + W{t), (1) 

which consists of two distinct parts. The static part 



Ho 



j=i veA 

3 3. 
i=l j = lj^i ^ reA 



reB 



(2) 



is simply the tight-binding Hamiltonian in the honey- 
comb lattice, where a J and (br) are, respectively, 
fermionic creation and annihilation operators on the lat- 
tice site r in the A {B) sublattice. The three vectors 

di = dcx, d2 = ^{-^x + VSCy), ds = ^(-e^ - VSCy), 

(3) 

connect an A-lattice site with its three nearest-neighbor 
(nn) 5-lattice sites and are given in terms of the distance 
d = Stt/Sa/S/c between nn sites, where k is the laser wave 
number (see Fig.jl]). Here, 7,7' > characterize the 
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FIG. 1: (Color online) Laser configuration to create the hon- 
eycomb lattice, which consists of two triangular sublattices 
(A, red dots, and blue dots). The vectors di, d2, and ds 
connect a site on the A sublattice to its nearest neighbors on 
the B sublattice. 

The time-dependent part of the Hamiltonian ([T]), 



W{t) = mn^ cos{nt) ( ^ r • palar + r • pblbr 



KreA 



reB 



(4) 

describes the harmonic shaking of the lattice in the direc- 
tion p with a driving frequency ft in the co-moving frame 
of Ref. [16]. As a consequence of the transformation to 
the co-moving frame, W{t) describes atoms of mass m 
experiencing a position-dependent sinusoidal force. 



B. Effective Hamiltonian 

The unavoidable complexity that arises when dealing 
with a quantum many-body system out of equilibrium 
has recently motivated the development of new theoret- 
ical tools, for example time-dependent density matrix 
renormalization group [17 , time-dependent dynamical 
mean field theory [18], and exact diagonalization [19^. 

However, for a periodically driven quantum system, 
the Floquet theory offers a simplified description of 
the system, in the form of a time-independent effective 
Hamiltonian, if the period T = 27r/r^ is the shortest time 
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scale in the problem [20]. In this limit, the atoms cannot 
follow the shaking motion adiabatically and remain thus 
at their average lattice position, albeit with renormalized 
hopping parameters. The system is thus considered to be 
in a stationary state and the knowledge of equilibrium 
physics can be employed. 

Let us consider the Floquet Hamiltonian defined by 
Hf = H{t) — ihdt^ where H{t + T) = H{t) is periodic in 
time [2^. The eigenvalue equation is then given by 

HF\^{q,t)) = e^\^{q,t)), (5) 

where is the quasienergy defined uniquely up to a mul- 
tiple of HQ. Any solution | (j){q^ t)) is part of a set of solu- 
tions ex.p{inQt)\(j){q,t)) with integer n, which all corre- 
spond to the same physical solution. Hence, the spectrum 
of the Floquet Hamiltonian possesses a Brillouin-zone- 
like structure ^20|. The interest therefore lies with the 



states in the first Brillouin zone, i.e. states with quasi- 
energies —hVt/2 < < 

The space in which the states \4>{q^t)) are defined is 
the composite of the Hilbert space spanned by square 
integrable functions on configuration space, \a{q))^ and 
the space of T-periodic functions. The state | (j){q^ t)) may 
be written down in an orthonormal basis in the composite 
space according to 

oo 

\<t>{Q,t)) = ^^c„,aexp[-iF(i) +mm]|a(9)), (6) 

n=0 a 

where Cn^a are coefficients to normalize | <p{q^t)) and the 
operator F{t) can be any T-periodic Hermitian operator. 
Therefore, we can conveniently choose F{t) to be F{t) = 
h-^ /q dtW{t'), such that H{t) - hdtF(t) = Hq. 



If the condition 

{a\q)\{exp[iF{t)] exp[i{n - n')nt] Hoexp[-iF{t)])T\ a{q)) <C HQ (7) 
is satisfied for any two states | a{q)) and | a\q)), then the eigenvalues are approximately 

= {^{q,t)\{HF)T\^{q,t)) - ^ co,a'Co,a{^\q)\ {exp[iF{t)] Ho exp[-iF{t)])T \ a{q)). (8) 

a, a' 



Here, {0{t))T = T'^ /q dt 0{t) denotes the time average 
of the operator 0{t) over the period T. The condition 
^ will hold for n ^ n' ^ if i^o is nearly constant during 
the period T, which is small if ft is large. In this case, 
states with different n do not mix. If Q is large enough, 
such that the condition (7|) also holds for n = n\ then the 
energy spectrum will split up into energy bands labelled 
by an index n, where the details within the energy band 
are determined by Hq. Because the states with a different 
index n are separated by an energy which is a multiple of 
h^} and because the spectrum possesses a Brillouin-zone- 
like structure, only the terms with n = need to be taken 
into account. The effective Hamiltonian i^eff, which gives 
rise to the same spectrum as the Floquet Hamiltonian, is 
then defined by [21] 

i^eff = (^exp[iF{t)]Ho exp[-iF{t)]^ 

CO \ 

^5-[F(t),i7o]n) . (9) 

' n=0 ^' I T 



Here, [F, G]^ denotes the multiple commutator, which is 
defined by [F, G]n+i = [F, [F, G]n] and [F, G]o = G. 

Effective Hamiltonians corresponding to Eq. ([9| have 
been derived for linear shaking of a one-dimensional lat- 
tice [22] and for elliptical shaking of a triangular lattice 
[23] . For the shaken honeycomb lattice studied here, the 



condition ^ is satisfied if 7 <C hTt^ and the effective 
Hamiltonian becomes 



3 3 . 

- XI ^Ij f «r+d,-d, + ^I^r+d,-d 

i=l j = l,jyi ^rGA veB 

-/if ^ajar + ^hlhA, 



veA 



(10) 



where the renormalized nn hopping parameters jj are 
given by 



7j = 7^0 



(11) 



and the renormalized nnn hopping parameters are given 

by 



7L- = I'Jo 



(di-dj) - p 



IT 



(12) 



(see Appendix [A] for detailed calculations). In these ex- 
pressions, Jo (x) denotes the zeroth order Bessel func- 
tion of the first kind, which shows a damped oscillation 
around zero. 
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In terms of the renormalized nn and nnn hopping pa- 
rameters, the diaginahzation of the effective Hamiltonian 
( 10 ) yields the dispersion relation 



eA(q) = Mq) + A|/(q)|, 



(13) 



where A = ± is the band index, and we have defined the 
functions 

/(q) = ^7jexp(-iq-d,) (14) 



and 



III. 



/i(q) = 2^7;,^.cos[q-(d,-d,)]. (15) 

i<j 



MERGING AND ALIGNMENT OF DIRAC 
POINTS 



In this section, the honeycomb lattice with anisotropic 
hopping is studied. In the first two subsections, only nn 
hopping is considered for illustration reasons. Indeed, 
this allows for a simple understanding of the main con- 
sequences of shaking on Dirac-point motion and dimen- 
sional crossover. In Subsec. |III C| we discuss how the 
picture evolves when nnn hopping is included, and the 
sign of the hopping parameters is investigated in Subsec. 



IIID[ Since the system with two nn equal hopping param- 
eters and a single independent one captures the essential 
features of the systems with three independent nn hop- 
ping parameters, we will focus on this system. The num- 
bering of the 7jS is chosen such that I72I = I73I = 72,3, 
which can be achieved by shaking in a direction parallel 
or perpendicular to di. 

Although the atoms in the optical lattice are charge- 
neutral objects, we shall adopt the language from 
condensed-matter physics and call a zero-gap phase with 
a pair of Dirac cones and a vanishing density of states at 
the band-contact points a semimetal, whereas a gapped 
phase is called band insulator. Furthermore, nnn hop- 
ping induces a metallic phase for small values of be- 
cause of an overlap between the two bands that yields a 
non- vanishing density of states at the energy level of the 
band-contact points. 



A. Merging of Dirac points 



zone, start to move in the g^y-direction along the verti- 
cal edges of the latter. This motion is depicted by the 
arrows in Fig.jsj^b). Even if the two Dirac points are 
no longer located at the high-symmetry points K and 
K\ they remain related by time-reversal symmetry, such 
that their Berry phases tt and — tt are opposite. This 
non-zero Berry phase topologically protects each of the 
Dirac points and thus the semimetallic phase remains 
robust until 72,3 = 7i/2, where the two points merge at 
a time-reversal invariant momentum, i.e. half of a re- 
ciprocal lattice vector [6]. In the present example, this 
point is situated at the center of the vertical edges of 
the first Brillouin zone, and the band dispersion becomes 
parabolic in the ^/-direction while remaining linear in the 
x-direction [see Fig.|2|b)]. The merged Dirac points are 
no longer topologically protected due to the annihilation 
of the opposite Berry phases. Consequently, a further in- 
crease of the shaking amplitude, which results in a further 
decrease of 72,3, leads to the opening of a gap between 
the two bands. Thus, the system undergoes a topological 
phase transition from a semimetal to a band insulator. 
This merging transition was also studied in a static setup 
in Ref. where the hopping amplitudes 7 were pro- 

posed to be modified by a change in the intensity of one 
of the lasers used to create the optical lattice. In contrast 
to this static setup, shaking the honeycomb lattice allows 
one to completely annihilate some of the nn hopping pa- 
rameters and to even change their sign. This sign change 
occurs at the zeros of the Bessel function [see Eq. (11)]. 
For an example system of ^^K atoms in a lattice created 
by lasers with a wavelength of 830 nm, which is shaken in 
the direction perpendicular to di, the situation 72,3 = 
is encountered for 



p = 180nm; n/27T = 6kHz, 



(16) 



which corresponds to the first zero of the Bessel func- 
tion. At this particular point, and if 7' = in addition, 
the system consists of a set of effectively decoupled hor- 
izontal bonds along which the atoms are solely allowed 
to hop. This yields two fiat bands at ±71 [see Fig.[2|c)] 
that may be viewed as the extreme limit of the band- 
insulating phase. Alternatively, one may view this situ- 
ation upon decreasing the value of 72,3 as a dimensional 
crossover from a 2D band insulator to a zero-dimensional 
(OD) system. A small non-zero value of 72,3 simply pro- 
vides a weak dispersion of these decoupled bands (not 
shown) . 



If the latice is shaken in the direction perpendicular to 
di (direction 1 in Fig.jl]), 71 remains equal to 7, whereas 
72 and 73 are renormalised to a smaller value. An in- 
crease in the shaking amplitude results in a decrease in 
72,3 =72 =73, which is depicted by the arrow M in 
Fig.jsja). When the hopping parameters change accord- 
ing to this arrow M, the energy spectrum evolves from 
Fig.|2|a) to Fig.|2|b). The Dirac points, originally sit- 
uated at the corners K and K' of the first Brillouin 



B. Alignment of Dirac points 

Another dimensional crossover, from 2D to ID, may be 
obtained if the lattice is shaken in the direction parallel 
to one of the nn vectors (direction 2 in Fig.[T]). Here, 
we choose di to maintain the symmetry 72,3 = 72 = 
73. In this case, both 71 and 72,3 are renormalized by 
Bessel functions, albeit with different arguments. Since 
all hopping parameters are renormalized, the trajectory 
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FIG. 2: (Color online) Energy dispersion for the shaken hon- 
eycomb optical lattice, with /c = 1, 7 = 1, and 7^ = 0. The 
labels qx and qy represent the x and y components of the 
momentum, respectively. The x and y axes have been chosen 
such that the nn vectors are given by Eq. ([3|. (a) The isotropic 
case, where 71 = 72,3- (b) The merged Dirac points, where 
72,3 = 71/2. (c) The zero dimensional case, where 72,3 = 0. 
(d) The aligned Dirac points, where 71 = 0. 
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for the same system of ^^K atoms mentioned above. Here, 
for illustrative purposes, a simplified trajectory of the 
system is depicted in Fig. [2] by arrow A, which corre- 
sponds to the motion of the Dirac points in reciprocal 
space as shown in Fig.jsjc). As 71 approaches zero, 
the Dirac points align in lines parallel to the x-axis at 
Qy = :^7^/^/3d and the energy barriers between the align- 
ing points are lowered. Consequently, when 71 = 0, 
the energy spectrum contains lines where the two energy 
band meet and the dispersion is linear, as is shown in 
Fig.[2|d). The dispersion relation ([T3| then reads simply 



eA(q) = 2A72,3 



cos I —qyd 



(18) 



and one clearly sees the ID character. Indeed, there is 
no dispersion in the g^a^- direct ion, as is also evident from 
Fig.[2|d), and the system may be viewed as completely 
decoupled ID chains in which the zig-zag arrangement 
is of no importance. In this particular limit, the sites A 
and B are therefore no longer inequivalent such that the 
unit cell is effectively divided by two, and the size of the 
first Brillouin zone is consequently doubled. The aligned 
Dirac points may thus, alternatively, be viewed as due to 
an artificial folding of the second (outer) half of the first 
Brillouin zone into its inner half. However, this aspect 
is very particular in that the Brillouin zone immediately 
retrieves its original size when 71 is small, but non-zero, 
or if nnn hoppings are taken into account. In both cases, 
one needs to distinguish the two different sublattices and 
one obtains a dispersion in the g^^^- direct ion. 

The actual behavior of the system for an increasing 
shaking amplitude is discussed in Sec.|V[ This behavior 
is more complicated because all three nn hopping pa- 
rameters are renormalized, which, beyond the alignment, 
leads to the merging of Dirac points and the opening of 
a gap also in the case of shaking parallel to one of the nn 
vectors. In the absence of nnn hopping, the OD limit can 
be reached in addition. 



FIG. 3: (Color online) (a) Phase diagram, showing the phase 
transition between the zero-gap semimetallic phase and the 
insulating phase, which happens at 71 = 272,3- Here, we have 
chosen 7^ = 0. (b) Dirac-point motion in the first Brillouin 
zone for a shaking direction perpendicular to di [direction M 
in the phase diagram (a)], (c) Dirac-point motion in the first 
Brillouin zone for a shaking direction parallel to di [direction 
A in the phase diagram (a)]. The contour plots depict the dis- 
persion of the isotropic system with an arbitrary color scale. 
The area with higher contrast is the first Brillouin zone. 



of the system in the phase space upon increasing the 
shaking amplitude is not a straight line, as was the case 
for shaking perpendicular to a hopping direction, and 
has a new feature: the alignment of Dirac points, which 
occurs for 71 = 0. The first zero of 71 is found at 



C. Next-nearest-neighbor hopping 

The major consequence of nnn hopping is to break 
particle- hole symmetry, as may be seen from Eq. (flsl). 



92nm; n/27T = 6kHz, 



(17) 



where non-zero values of yield ex{q) 7^ — e_A(q). Its 
relevance depends sensitively on the shaking direction, 
because of the different renormalization of the nn hop- 
ping parameters. The band structure with nnn hopping 
included is depicted in Fig. [4] for different shaking direc- 
tions. 



1. Shaking in the direction perpendicular to di 

In the case of a shaking perpendicular to di, only 72,3 
are decreased, whereas 71 = 7 remains the leading en- 
ergy scale in the band structure ^32^. The band struc- 
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ture for the unshaken lattice is depicted in Fig.|4ja) for 
7Y7 = 0.1, and one notices that the main features of 
the band structure, namely the Dirac points, are unal- 
tered with respect to the case 7' = in Fig.[2|a), apart 
from the flattening of the upper band as compared to 
the lower one. When approaching the merging transi- 
tion 72,3 = 7i/2, the value of which is determined by the 
zeros of /(q) in Eq. (14) and that therefore does not 



depend on the nnn hopping parameters, the band width 
remains dominated by the largest hopping parameter 71 , 
such that the band structure [Fig.[4]^b)] at the transi- 
tion is essentially the same as in Fig.^b) for 7' = 0. 
In the OD limit, with 72,3 = the originally flat bands 
[Fig.|2|c)] acquire the weak dispersion of a triangular lat- 
tice as a consequence of the non-zero nnn hopping param- 
eters. However, as expected from the above arguments, 
the dispersion is on the order of 7' and thus small as 
compared to the energy separation ~ 271 = 27 between 
the two bands. 
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2. Shaking in the direction parallel to di 



In contrast to a shaking direction perpendicular to di, 
nnn hopping has more drastic consequences if the lat- 
tice is shaken in the direction parallel to di. In this 
case, all nn hopping parameters are decreased, and the 
relative importance of nnn hopping is enhanced. No- 
tice further that the nnn lattice vectors ±(d2 — d3) are 
now perpendicular to the shaking direction such that 

3 — 1^ — O-^T remains unrenormalized. Also in this 
case, the system is approaching the ID limit, with 71 = 
[see Fig.|4]^d)]. However, in contrast to Fig.[2|d), the 
chains remain coupled by nnn hopping that yields a dis- 
persion in the g^j-direction. Furthermore, as mentioned 
above, the A and B sites are now not equivalent from a 
crystallographic point of view, such that the outer parts 
of the first Brillouin zone cannot be folded back into the 
inner one, as may be seen from Fig.pF 



Finally, for particular values of the shaking ampli- 
tude in the direction parallel to di, the nn hopping 
parameters can be decreased in such a manner as to 
render 72,3 more relevant. In this case, the two bands 
can overlap in energy, as depicted in Figs.[4]^e) and 
for p = b.2{h/mVtd)ex (in which case 71 ^ 72,3) and 
p = A.S{h/mVtd)ex (with 72,3 ~ 0), respectively. In the 
latter example there are no band contact points, in spite 
of the overlap between the two bands, and the system 
would be in an insulating phase if nnn hopping terms 
were not taken into account. This overlap in energy be- 
tween the two bands yields a non-zero density of states 
at any energy, such that the semi-metallic (or insulat- 
ing) phase vanishes and yields, at half-filling, a metallic 
phase with particle and anti-particle pockets in the first 
Brillouin zone. 



FIG. 4: (Color online) Energy dispersion for the shaken hon- 
eycomb optical lattice, with A; = 1, 7 = 1, and 7^ = 0.1. 
The labels qx and qy represent the x and y components of 
the momentum, respectively. The x and y axes have been 
chosen such that the nn vectors are given by Eq. (|3|. (a) 
The homogeneous case, where 71 = 72,3- (b) The merged 
Dirac points, where 72,3 = 7i/2. (c) The zero dimen- 
sional case, where 72,3 = 0. (d) The aligned Dirac points, 
where 71 = 0. (e) An example of the metallic phase with 
p = 5.2{h/mQd)ex. (f) Another example of the metallic phase 
with p = 4:.8{h/mQd)ex. 



D. The signs of the hopping parameters 



As already alluded to in the previous sections, the 
shaking of a honeycomb lattice can lead to a sign change 
of the hopping parameters. Quite generally. Fig. |5] shows 
that changing the relative signs of the nn hopping pa- 
rameters results in a translation of the energy spectrum 
in momentum space. This effect was also mentioned in 
Ref. [4 . Indeed, the relative signs determine at which of 
the four time-reversal invariant momenta in the first Bril- 
louin zone the merging of Dirac points and the semimetal- 
insulator transition take place when |7i| = 2|72,3|. How- 
ever, the sign change of the nn hopping parameters can 
be transformed away by a gauge transformation . Nev- 
ertheless, the sign of the nnn hopping parameter is im- 
portant, since it determines whether the upper or the 
lower band is flattened. 
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FIG. 5: (Color online) Phase diagram and contour plots of 
the energy bands, showing the effects of the renormalized nn 
hopping parameters 7^ . (a) to (h) Contour plots of the en- 
ergy bands, where the color scaling is arbitrary and the first 
Brillouin zone is the area with higher contrast. The value of 
the nn hopping parameters for each contour plot are given by 
the position of the corresponding letter in the phase diagram 
and 7^ = 0. The dark regions indicate energies close to zero, 
whereas brighter regions are further away in energy from the 
Fermi level at half filling. 



IV. INTERACTIONS 

Until now, we have considered single-component 
fermionic atoms and, due to the Pauli principle, the ab- 
sence of s-wave interaction naturally results in an ideal 
Fermi lattice gas, albeit with an unusual band structure. 
By trapping two hyperfine states of the fermionic atoms, 
Hubbard-like interaction terms arise. 

Hint = X] X] ^ 4,a4,a'^r,a'^^^r,a 



r^A a, a' 



u 



(19) 



where the fermionic operators now acquire an additional 
spin index a = {t, i}, which is summed over, and U is 
the interaction energy. Naturally, in order to be able to 
apply the Floquet theory in the presence of interactions, 
the Hamiltonian Hq in Eq. ^ must now be replaced by 
H = Ho-\- Hint' This is the case in our study because we 
investigate the weak-coupling limit with /7 <C 7. Since 
the interaction term commutes with the shaking, it is not 
renormalized, similarly to the on-site energy term propor- 
tional to /i in Eq. ([2|. However, it has been shown that 
complications may arise when a multiple of the energy 
U is in resonance with a harmonic of hQ^ mhQ = n/7, 
for integer m and n. Whereas the limit [25 m ^ n is 
not considered here because it is in contradiction with the 



small- large-frequency limit, critical resonances may oc- 
cur for m n [26j. Nevertheless, it has been shown in 
Ref. [26] that these resonances, which occur in higher- 
order perturbation theory, are strongly suppressed in the 
large-m limit. 

In the weakly-interacting regime considered here, the 
ground state is adiabatically connected to that of the 
non-interacting system, with no broken symmetry. First, 
we use the Fourier transform of the creation and annihi- 
lation operators, a^^a = exp(iq-r)aq^cr, to find 
the Hamiltonian in momentum space. Within a Hartree- 
Fock theory, we introduce a mean-field decoupling of the 
interaction terms. 



(20) 



(^q2,cr'^q3,cr')<^ql,cr<^q4,cr " (<^q2,cr' ^q4,cr )<^ql,cr^q3,cr' 
+ ^q2,cr'^q3,cr' (<^ql,cr^q4,cr) " <^q2,cr' ^q4,cr (<^ql ,cr^q3,cr' ) 
-(42,a'«q3,a')(«ql,a«q4,a) + (^^q2,a' ^q4,a ) (aql,aaq3,a' ) ^ 

such that the expectation values of both sides are equal. 
For the mean value we take 

(^q,a^q',a') = (^q,aV,^') = ATriq^^^q^q/ ^^,^/ , (21) 

where M is the number of sites per sublattice, nq^cr is 
the density of atoms with momentum q and spin index 
cr, and ^a.ot' is the Kronecker delta. We then obtain the 
mean-field Hamiltonian 



MF 



H, 



UMn^ 



eff ■ 



Un 

cr,q 



(^q,cr^q,o 



^q,a^q,a). 



where the total density is defined by n = ^ ^ \ 



(22) 



The Hamiltonian (22) may be rewritten in a matrix 
form: 



Hmf 



UAfr. 



(23) 



cr,q 



^(M,q) /(q) 

/*(q) /^(M,q) 



^q,cr 



where we have introduced the functions 

3 3 



/i(/i,q) 



Un 



-/i-y^ ^ exp[-zq-(d,-d,)], (24) 

i=l j = ljy^i 



and /(q) is defined in Eq. (14) The Hamiltonian (23) 
can then be diagonalized by the unitary operator 



which yields 



V2 



i/(q)/l/(q)l 

—i 



(25) 



8 



H 



MF 



8~~ 



MM,q)-|/(q)l 

/l(Ai,q) + |/(q) 



Because the c and d quasiparticles are free, the partition function corresponding to the Hamiltonian (26) reads 
Z = exp (^log |l + exp [ - p{h{^i,ci) - |/(q)|)] | + log |l + exp [ - /?(/i(/x,q) + |/(q)|)] 



cr,q 



(26) 



(27) 



where {3 = {kBT) ^ with denoting the Boltzmann's constant and T the temperature. 
The total number of particles N is given by (l//3)91og and one obtains 



cr,q ^ 



+ 



exp [/?(/i(/x,q) - |/(q)|)] 1 + exp [p{h{^i,a,) + |/(q)|)] 



(28) 



Since the expression inside the sum does not depend on spin, summing over a yields a factor 2. One recognizes in 
Eq. (28) the Fermi-Dirac distribution function Nfd{x) = [1 + exp(x)]~^. The number of particles N is related to the 
density n, which is defined here as the number of particles per lattice site, i.e. n = N/2M. Converting the sum over 
q into an integral, the following self-consistent equation for the density is derived 



n(/i) 



IBZ 



(29) 



where the integral is restricted to the first Brillouin zone, the surface of which is Vibz- 



In Fig.[6]^a), the density n(/i) is plotted for several val- 
ues of 72,3/71. For the isotropic case, 72,3/71 = 1, the 
result of Zhu et al. is reproduced [27 . For the OD limit, 
the fiat line due to the gap in the spectrum is clearly vis- 
ible at the chosen temperature. Fig.[6](b) confirms that 
repulsive interactions lead to a lower density than in a 
system without interactions for the same chemical po- 
tential. Fig.[6jc) agrees with the observation that the 
nnn hopping breaks the particle-hole symmetry. This ef- 
fect is also visible in Fig.[6]^d), where the dependence of 
the density on the chemical potential is calculated for a 
shaking vector where the system is in the zero-gapped 
semimetallic phase for 7' = and in the metallic phase 
for y = 0.1. 



V. POSSIBILITIES FOR EXPERIMENTAL 
OBSERVATION 



Honeycomb optical lattices have recently been realized 
experimentally, although the existing setups have only 
been used to investigate bosonic atoms [8l|28]. Shaking of 
a lattice has been experimentally implemented in a one- 
dimensional one by a periodic modulation on the position 
of the refiecting mirrors [14 . For a honeycomb lattice, 
the shaking could be realized by means of an acousto- 
optical device, as proposed for a triangular lattice in Ref. 
B3j. 

The magnitude of the nn hopping parameter 7 in a 



honeycomb optical lattice has been evaluated in Ref. 



7 ^ 1.861^;^? 



Vr 



ErJ 



3/4 



exp -1.582 




(30) 



in terms of the recoil energy Eji = h^k'^ /2m and the mag- 
nitude of the potential barrier between nearest-neighbor 
lattice sites Vq. The magnitude of the nnn hopping pa- 
rameter 7' is not yet known, but could be determined 
from numerical band structure calculations. In a typical 
experimental situation, we expect the ratio 7V7 to be 
in the 5 — 10% range, in agreement with the parameter 
chosen in the discussion of Sec. lIII Cl 

In a typical experiment, the shaking amplitude would 
be increased from zero to a finite value. Fig.[7| shows in 
which order the system goes through the different phases 
and dimensions upon increasing the shaking amplitude. 
Here, also the values of the shaking amplitude required 
for the dimensional crossovers are given for the same sys- 
tem as discussed in Sec. HI C and 7' = O.I7. If the shak- 
ing direction is perpendicular to one of the nn vectors, 
the system will be in the gapped insulating phase beyond 
a certain value of the shaking amplitude, since the Bessel 
function crosses the value 0.5 only once and never obtains 
the value -0.5. If the shaking is parallel to one of the nn 
vectors, the system will be in the metallic phase beyond 
a certain value of the shaking amplitude. Nevertheless, 
it is still possible to induce a merging of Dirac points 
and to open up a gap in the spectrum, since the nn hop- 
ping parameters are renormalized such that the value of 
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FIG. 6: (Color online) Density n as a function of the chemical 
potential /x. Unless specified otherwise in the figure, the nn 
hopping parameters 72,3 = 7i = 7 = 1, the nnn hopping 
parameter 7^ = 0, the interaction strength [7 = 0, and the 
inverse temperature /3 = 20. (a) Effect of the renormalization 
of the nn hopping parameters, (b) Effect of the interaction 
strength U for the isotropic case, (c) Effect of the nnn hopping 
parameter 7^ = 0.1 in the shaken lattice, (d) The metallic 
phase. For the isotropic cases, p — 5.2(h/mQd)ex, whereas 
for the OD cases p = 4.8{h/mQd)ex. These systems are in the 
metallic phase for j = 0.1, whereas for 7^ = they are in the 
semi-metallic and the insulating phase, respectively. 



one of them will in general differ from that of the other 
two. However, whether the system is actually driven into 
an insulating phase or remains metallic depends on the 
precise value of the ratio 7V7- 
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FIG. 7: (Color online) Overview of the different phases as a 
function of the shaking amplitude. The bottom scale gives 
the size of the argument of the Bessel function, whereas the 
values for p in nm correspond to the optical lattice discussed 



in Sec. [mCl with 7' = O.I7, Q/27v = 6kHz, k = 830nm, and 
containing °K atoms, (a) Shaking perpendicular to one of 
the nn hopping directions, (b) Shaking parallel to one of the 
nn hopping directions. 



In experiments, an overall harmonic trapping potential 



Vtrap(r) 



(31) 



where entrap is the trapping frequency, and r is the posi- 
tion measured from the center of the trap. By apply- 
ing the local density approximation (LDA), one finds 
that the chemical potential evolves radially according to 

p^ p - Vtrap(r^). 

Fig.jsj^a) shows the density profile for several ratios of 
72,3/71, without nnn hopping or interactions. The case 
with 72,3/71 = 0, when the system is in the extreme limit 
of the band insulating phase, can be well distinguished 
from the other cases. Fig.jsjb) shows that stronger in- 
teractions lead to a higher density away from the cen- 
ter of the trap. This effect becomes visible when the 
density starts to deviate from one particle per lattice 
site. Next-nearest-neighbor hopping leads to a higher 
density at the edge of the cloud compared to the case 
without nnn hopping, which can be seen from comparing 
Figs. [s]^ a) and (c) and from Fig.jsjd). The latter shows 
the effect of nnn hopping on the density profile for the 
case where the nnn hopping gives rise to the metallic 
phase for two different shaking vectors. In the first case, 
p = b.2{h/mQd)ex^ which gives 71 ^ 72,3, such that 
without nnn hopping, the system is in the zero-gapped 
semi-metallic phase and the Dirac points are located very 
close to the corners of the first Brillouin zone. In the sec- 
ond case, p = A.8{h/mVtd)ex^ which results in 72,3 ~ 0, 
such that without nnn hopping the system is in the in- 
sulating phase and the two energy bands are almost flat. 

We emphasize that, in the present paper, we only dis- 
cuss weak correlations that adiabatically affect the den- 
sity. However, when increasing further the onsite inter- 
action, one may expect correlated phases with inhomo- 
geneous density, even at half-filling. A detailed study of 
these correlated phases is a vast research issue that is yet 
ongoing and that is beyond the scope of the present pa- 
per. Here, we only provide a glimpse on how the density, 
which we discussed above in the weak-coupling limit, may 
evolve in view of some phases studied in the literature. 
Mean- field calculations indicate a transition to an antifer- 
romagnetic state above a value oiU 0^2.2 [9 , whereas 
more sophisticated quantum Monte- Carlo calculations 
indicate an intermediate spin-liquid phase between the 
semimetal and the anti-ferromagnetic phase [10]. The 
spin-liquid phase may be viewed as a Mott insulator with 
a charge localization on the lattice sites, and recent slave- 
rotor calculations indicate that such spin-liquid phases 
dominate the phase diagram for 71 > 72,3 [11 , which is 
the parameter range where the Dirac points would merge 
in the absence of interactions. The precise transition 
between the weakly-interacting liquid phases and these 
strongly-correlated Mott insulators could in principle be 
determined with the help of the above-mentioned density 
measurements. 

A more promising technique for detecting Dirac-point 
motion is momentum-resolved Raman spectroscopy. This 
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FIG. 8: (Color online) Density n as a function of the dis- 
tance from the trap's centre r = |r|, which is expressed in 
units of the nearest-neighbor distance d. The trapping fre- 
quency has been chosen such that the trapping potential is 
given by Vtrap(r) = 0.00l7r^/(i^. The chemical potential \i 
for each case has been chosen such that the density at the 
trap's centre n is one particle per site. This corresponds to 
half-filling, since we consider 2 species of fermions. Unless 
specified otherwise in the figure, the nn hopping parameters 
72,3 = 7i = 7 = 1, the nnn hopping parameter ^' — 0, 
the interaction strength [7 = 0, and the inverse temperature 
/3 = 20. (a) Effect of the renormalization of the nn hop- 
ping parameters, (b) Effect of the interaction strength JJ for 
the isotropic case, (c) Effect of the nnn hopping parameter 
7^ = 0.1 in the shaken lattice, (d) The metallic phase. For the 
isotropic cases, p — b.2{h/mQ.d)ex^ whereas for the OD cases 
p = A.S(h/mQ.d)ex. These systems are in the metallic phase 
for 7^ = 0.1, whereas for 7^ = they are in the semi-metallic 
and the insulating phase, respectively. 



technique has been proposed as an equivalent of angle- 
resolved photoemission spectroscopy for cold atom sys- 
tems [2F. It has not yet been realized experimentally, 
while momentum-resolved radio-frequency spectroscopy, 
which is a very similar technique, has already been im- 
plemented [30 . Notice further that another very sim- 
ilar technique, momentum-resolved Bragg spectroscopy, 
has been applied to ultracold bosonic atoms in a static 
optical lattice by Ernst et al. [31]. Momentum-resolved 
spectroscopy can allow us to indirectly visualize the band 
structure. In momentum-resolved Raman spectroscopy, 
two laser pulses with frequencies ui and 002 are irradi- 
ated upon the system. If the frequency difference is in 
resonance with a transition uj^f between atomic hyper- 
fine states, uji — UJ2 = ^hf^ some atoms are excited in a 
second-order process to the higher hyperfine state. Then, 
with state-selective time-of-flight measurements, the dis- 
persion of the atoms in the new state are measured, from 
which the dispersion of the original atoms can be derived. 
When the atoms are confined in a trapping potential and 
the laser pulses are focussed on the center of the trap, the 
quality of the results obtained by Raman spectroscopy is 
comparable to those of a homogeneous system [W. Fur- 



thermore, Raman spectroscopy yields better results for 
a system with strong interactions compared to standard 
time-of-flight measurements [29 . 

Notice that momentum resolved Raman spectroscopy 
was originally proposed to be applied to a gas of ultra- 
cold fermionic atoms at equilibrium and not for a shaken 
lattice. We therefore discuss, in this final paragraph, why 
we think that this technique may also be applied to the 
present case. Naturally, as long as the frequencies of the 
additional lasers in the Raman-spectroscopy setup are 
small with respect to the shaking frequency, uji^uj2 ^ 
even the full system satisfies the condition ([7| for the va- 
lidity of Floquet theory. As in the case of interactions, 
one needs, however, to avoid resonances between the dif- 
ferent laser frequencies that could become critical [2F. 
The opposite limit, in which the laser frequencies and 
that of the hyperfine transition are larger than the shak- 
ing frequency, is more delicate. However, even then, the 
shaken system remains at quasi-equilibrium as long as 
the intensities of the lasers used in Raman spectroscopy 
are weak, such that they only constitute a small per- 
turbation. The atomic dynamics probed even at high 
frequencies is therefore still that of the atoms at quasi- 
equilibrium, with the band strucure obtained from Flo- 
quet theory. Furthermore, in the experimental studies 
by Zenesini et al. [14], time-of-flight measurements were 
used to determine the momentum distribution of bosonic 
atoms in a shaken lattice. Apart from the time-scale con- 
siderations, there are also some length scales that need 
to be taken into account. There are indeed two require- 
ments for the correct size of the focus of the laser beams. 
On the one hand, it needs to be larger than the lattice 
spacing, such that sufficiently many atoms can be ex- 
cited, while on the other hand the focus of the beams 
should be small enough, in order to have an approxi- 
mately flat trapping potential inside the focus area. In 
addition, choosing the length of the pulses could possibly 
be a problem, since for shorter pulses, the excited atoms 
will be less affected by the lattice potential, whereas for 
longer pulses more atoms can be excited, leading to a 
stronger signal. 



VI. CONCLUSIONS 

In conclusion, we have investigated the band engineer- 
ing of fermionic atoms in an optical honeycomb lattice 
with the help of a periodic shaking of the lattice. If the 
shaking frequency Q is large enough, i.e. if HQ constitutes 
the largest energy scale in the system, the Floquet theory 
may be applied and the system is at quasi-equilibrium in 
the sense that the atoms cannot follow the rapid motion 
associated with the shaking. Depending on the direction 
of the shaking, one may render the hopping amplitudes 
in the quasi-static lattice anisotropic, due to a renormal- 
ization of the nn and nnn hopping parameters by Bessel 
functions that go through zero and change sign. As a 
consequence, dimensional crossovers can be induced in 
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absence of the nnn hopping. For a shaking direction par- 
ahel to one of the nn vectors (such as e.g. di), one can 
make one of the nn hopping parameters vanish, 71 ^ 0. 
The system then undergoes a transition from 2D to ID, 
while the Dirac points ahgn simultaneously. Shaking in 
the perpendicular direction (± di) allows one to decrease 
two nn hopping amplitudes simultaneously while main- 
taining 71 unrenormalized. In this case, a dimensional 
crossover from 2D to OD is induced for 72,3 0, leading 
to two fiat energy bands, beyond the merging of Dirac 
points [11 [6], which occurs at |7i| = 2|72,3|. A nonzero 
value of 7' breaks the particle-hole symmetry and leads 
to a coupling among the ID chains and the OD dimers, for 
the 71 = and 72,3 = cases, respectively, and thus to a 
weak 2D dispersion. The merging and the alignment of 
Dirac points, however, are not affected. Moreover, for a 
shaking direction parallel to di, one pair of nnn hopping 
amplitudes [±(d2 — d3)] remains unrenormalized, and its 
relative importance is thus enhanced when compared to 
the decreasing nn hopping amplitudes. In this limit, be- 
yond the semi- metallic and the band-insulating phases, a 
novel metallic phase can appear that consists of particle 
and hole pockets with a non-vanishing density of states 
even at half- filling. 

Furthermore, we have investigated the role of weak 
repulsive on-site interactions. The resulting ground 
state is then adiabatically connected to that of the non- 
interacting system, and we have self-consistently calcu- 
lated the dependence of the atomic density on the (local) 
chemical potential. The density profiles of the different 
phases, e.g. the gapless semimetal or the gapped band 
insulator, and the different dimensionality may be mea- 
sured experimentally by in-situ density measurements. 
Moreover, momentum-resolved Raman spectroscopy may 
be a promising technique to measure the band structure 
associated with these different phases. 
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Appendix A: Effective Hamiltonian 

For the studied case, H{t) — hdtF{t) = Hq, where Hq 
was given in Eq. Since the nn hopping is usually 



larger than the nnn hopping, 7' < 7, and since we let 
the chemical potential be in the range —27 < fi < 27, 
the dominant energy scale in the Hamiltonian Hq is 7. 
Therefore, if 7 <C hVt^ the condition ^ is satisfied and 
the Floquet theory may be applied. 

In the general case, the effective Hamiltonian is given 

by m 



n=0 



(Al) 



where for the shaken honeycomb lattice, we choose 



F{t) 



sin(r^t) [ r • polar + r • pblbr 



\reA reB 

Using the (nonvanishing) commutation relations 

[al ttr', ^i+d,- ^r] = - ^i+d,- ^r' (^r',r, 
[al, a^'^al ar+d—d^] =^1, ar+d—d^ ^^r',r 

- al a^' ^r',r+di-dj, 



(A2) 



(A3) 



which are valid for both fermionic and bosonic creation 
and annihilation operators, and equivalent ones for the 
creation and annihilation operators on the B sublattice, 
the multiple commutator in Eq. (Al) becomes 



3 



sin(l]t) 



j=i veA 
3 3 

^'E E (E[(d,-di)-pra]:ar+d.-d, 

i=l j = l,jyi ^ reA 



^[(d,-d,).pr6t6r+d.-d,)| 



veA 



reB 



(A4) 



Finally, after performing the time average and evaluat- 
ing the sum over n, the effective Hamiltonian (10) is ob- 
tained. 
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